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I. INTRODUCTION 



Recently a program has been proposed to study neutrino masses This program is 
based on the following assumptions. 

Assumption 1: There are three right-handed Dirac neutrinos in addition to the three 
known left-handed neutrinos. 

Assumption 2: The lepton mass matrices have the same structure as the quark mass 
matrices, as described by Lehmann et al. @. 

In carrying out this program, one of the technically most difficult aspects is to treat the 
Mikheyev-Smirnov-Wolfenstein (MSW) effect for the sun [Q . It is the purpose of this paper 
to discuss this problem. 

The difficulty arises in the following way. On the basis of the Super-Kamiokande data on 
atmospheric neutrinos ||, the muon neutrino oscillates with another neutrino which is 
not the electron neutrino v e . From Assumption 1, this neutrino must be the tau neutrino v T . 
In other words, and v T are coupled significantly. Again from Assumption 1, the electron 
neutrinos v e from the sun must therefore oscillate with either or v T . Since and v T are 
coupled significantly, the electron neutrinos v e from the sun must oscillate with both and 
v T . It is accordingly necessary to understand the MSW effect with all three generations of 
neutrinos u e , and v T . 

While there is an extensive literature on the MSW effect for two generations of neutrinos 
HH, relatively little is known for the case of three generations @. One of the major 
differences between these two cases is the following. For two generations, the MSW effect 
occurs predominantly in the vicinity of a particular density in the sun. For this reason, 
the important parameter is the rate of change of the density at this particular density. In 
contrast, for three generations, the important contributions to the MSW effect can come 
from several regions. When the parameters, such as the neutrino masses, vary, these regions 
not only move but also merge and separate, features that are absent for the much simpler 
case of two generations. 



II. FORMULATION OF THE PROBLEM 



There are many different ways to formulate the problem of the MSW Q effect in the sun 
for three generations of neutrinos. These different formulations are of course closely related 
to each other; it has been found that, at least for the present program [I], the following one 
is most appropriate. 

In a current basis, the charged lepton currents are 



(2.1) 



where 



^3 



and 



denote the neutrino and charged lepton fields. The neutrino mass matrix M 
vacuum expectation value of the Higgs field in the neutrino-Higgs coupling, is 



from the 
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z/lMi/r + h.c. 



(2.2) 



Let p be the momentum of the neutrino, the value of which is taken to be much larger 
than the neutrino masses, which are the eigenvalues of M. Typically, p is in the range from 
0.4 MeV/c to 10 MeV/c, while the neutrino masses are believed to be no more than about 
1 eV; hence this assumption is well satisfied. Under this assumption, it is M 2 that enters in 
the MSW effect. Let 
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(2.3) 



be real and symmetric. Note that My are not defined, only M 2 j = (M 2 )ij. Of course the 
eigenvalues of this matrix M 2 are the squares of the neutrino masses. 

In terms of this M 2 , the MSW effect for neutrino oscillations is described by the coupled 
ordinary differential equation || 
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where 



D(r) = V2G F N e 



(2.4) 



(2.5) 



with Gp the Fermi weak-interaction constant and N e (r) the electron density at the point r. 

So far this formulation is general. In view of the observation in Sec. I that, for three 
generations of neutrinos, important contributions to the MSW effect can come from several 
values of r, we propose to choose a suitable function D(r). On the one hand, this D(r) 
should describe approximately the electron density in the sun; on the other hand, this D(r) 
should be sufficiently simple so that Eq. (|2.4| ) can be solved explicitly. This choice is made 
in the following way. 

First, since most neutrinos are produced by nuclear reactions Q near the center of the 
sun, the simplifying assumption is made that the neutrinos are produced at the center of 
the sun. This has the consequence that the r in Eqs. (|2.4| ) and fl2.5|) is the radial distance 
from the center of the sun. This is the reason for the notation r. 

Secondly, the density of the sun as a function of the radial distance is fairly accurately 
known [|n|, as shown by the solid curve in Fig. 1, where the vertical axis is the density in 



g/cm 3 in a logarithmic scale, and the horizontal axis is r/R Q . The radius of the sun is 

R Q = 6.96 x 10 8 m. 

From this solid curve in Fig. 1, it is seen that the solar density is approximately exponential, 
as shown by the dotted line of Fig. 1. This fit is reasonably good, to an accuracy of roughly 
15%, for 

0.05 < r/R Q < 0.9. 
In other words, D(r) is taken to be given by 

D(r) = D(0)e- r/ro . 



(2.6) 
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FIG. 1. Density profile of the sun [1C], with exponential approximation (dashed). 



From the dotted line of Fig. 1, the following values are obtained 
£>(0) = 0.0458 km" 1 , r ~ 0.1 x R Q . 
The problem is therefore to solve the MSW differential equation, for r > 0, 
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III. SOLUTION OF THE DIFFERENTIAL EQUATION 

The differential equation (|2.8|) can be solved explicitly in terms of generalized hyper- 
geometric functions. This is to be carried out in the present section. Indeed, this is the 
underlying reason for the choice of the exponential function (|2.6| ). 

Let 

u = r/r + u , (3.1) 

with 

M = - ln[D(0)r ] 8.07. (3.2) 



3 



Then 
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(3-3) 



Note that u = corresponds to 

r = —r u ~ 0.81i? Q . 
In view of the form of Eq. (|3.3|), it is convenient to diagonalize the 2x2 matrix 



(3.4) 



2p 



Ml Ml 
Ml Ml 



Let 8q be the angle of rotation that diagonalizes this 2x2 matrix, i.e., 
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with 



uj 2 < uj 3 . 



(3.5) 



(3.6) 



The reason for this Eq. (|3.6|) is that, if u 2 = 0U3, then the problem reduces to neutrino 
oscillations with only two neutrinos. Define also 



-Ml, 
2p 11 ' 



and 



ipi(u) = 4>i{u) 

then Eq. ( |3.3|) becomes 
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If X2 or is zero, then once again this reduces to neutrino oscillations with two genera- 
tions. By reversing the signs of ip 2 (u) and/or ip 3 {u) if necessary, it follows that, without loss 
of generality 



Xj >0, for j = 2, 3. 
Let fj,i, fi 2 and n 3 be the eigenvalues of the 3x3 matrix 



(3.10) 
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X2 X3 

X2 uj 2 ; 

then these three /i's are the squares of the neutrino masses multiplied by r /(2p). Note that 

3 3 

J2 UJ J = J2l x r ( 3 - n ) 

j=i j=i 

For cases of interest |], the fij will typically be large, fij ^> 1. This is of course not a 
question of units, but rather a reflection of the physical situation of having a large number 
of oscillations within the solar radius. 

The interlacing property implies that, by renumbering the /i's if necessary, 

< \i\ < uj 2 < /i2 < ^3 < A*3, (3.12) 

again assuming that there is no reduction to two-neutrino oscillation. 

The following properties are easily verified: With fixed u 2 , u> 3 , X2 and x 3 , 

(i) ^>0 for j = 1,2,3. 

(ii) As uji — > oo, 

Hi — > w 2 , /i 2 — > w 3 , and // 3 — > oo. 

(iii) As — > — oo, which is physically not possible, 

fit — > — oo, /i 2 — ► ^2, and // 3 — > cu 3 . 



A third-order ordinary differential equation can be obtained from Eq. ( |3.9| ) as follows. 
Define xpn, ip 22 and ^33 such that 

\ du J \ du J 
4>2 = X2 (i-^ -uj 3 ^j ip 22 , 

"03 = X3 (i^ ^2 J ^33- (3.13) 
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These ^ 22 and ^33 are of course not unique. It then follows from Eqs. ( |3.9| ) and ( |3.10| ) 
that 



(3.14) 



Therefore, 



^22 = fax + K 22 e~ l " 2U + K 23 e~ l ^ u , 
fa 3 = fax + K 32 e-^ U + K 33 e~^ u . 

In particular, 

^22 - (K 23 - K 33 )e~^ u = fa 3 - {K 32 - K 22 )e-^ u . 

Define this function of Eq. (|3.16|) to be ip, then it follows that 

fax = if; - K 22 e-^ u - K 33 e~^ u , 
fa2 = ^ + (K 23 - K 33 )e~^ u , 
fa 3 = iP + (K 32 - K 22 )e~^ u . 



(3.15) 



(3.16) 
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Substitution into Eq. ( |3.13| ) then expresses fa, fa and ip 3 in terms of the single function fa 
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fa = X3 (i j- - wa) ^. 



(3.18) 



Furthermore, this t/> is unique because of (|3.6|) 

To obtain the third-order ordinary differential equation for if), it only remains to substi- 
tute ( |3.18| ) into the first equation of (|3] 
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From the definition of fix, fi 2 and fi 3 , this can be rewritten as 
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To put this equation into the form of the differential equation for the generalized hypergeo- 
metric function, let 



z = ic 



(3.21) 



then the differential equation ( |3.20| ) for if) is 
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This is the differential equation for the generalized hypergeometric function 2 F 2 — see, for 
example, p. 184 of [IT|. Three linearly independent solutions of this third-order differential 
equation ( |3.22| ) are 
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where Kx, K 2 and i^ 3 are arbitrary non-zero constants. For the convenience of numerical 
calculations, we make the following choice: 



In 

Kx = [{ijJ 2 - fll){U)Z - fll)(fl2 - pii)(A*3 - Ml)] , 
if 2 = [(CJ 2 - fl2)(u 3 - fl 2 ){fll - fl 2 ){flz ~ fl 2 )\ 
if 3 = [(ufc - fl 3 ){0J 3 - /U 3 )(/il - //3j(A*2 - /is)] 

The function 2 F 2 can be defined by the following series expansion[] 0] 
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where the Pochhammer symbol is defined as 

(a)k = ol{ol + 1) • • • {ol + k — 1) 



(3.24) 



(3.25) 



(3.26) 



Since the Taylor series of 2 F 2 about z — 0, Eq. (|3.25|) , is always convergent, the differenti- 
ations specified in Eqs. ( |3.18| ) can be carried out term by term on these Taylor series. Thus 
these derivatives of can be easily expressed in terms of 2 F 2 with parameters slightly 
modified. The result is 



^3 = Ci4 1) + C 2 4 3 ) +C 3 4 3 3 \ 



(3.27) 



where 



From this form, it is also rather obvious how the confluent hypergeometric function 1F1 emerges 
as the solution when one of the a's equals one of the p's, i.e., for the case of two generations. 
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In Eqs. ( p.27| )) the constants Ci, C 2 and C 3 are determined by the boundary conditions 
Eqs. (|2.9| ). For the exponential matter distribution (|2.6|), this gives completely the solution 
of Eq. (|2.4j) , which describes the MSW effect in the sun. 

We recall that by Eq. e~ u = r D(0)e- r/r ° . Thus, the prefactors e~^ u are propor- 

tional to the exponentials describing propagation in vacuum, e 

-i«r/ro_ In fact 

as r — > 00, 

all these functions approach a simple limit: 
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While this solution ( |3.27| ) is specifically for three generations of neutrinos, generalization, 
if desired, to N generations is completely straightforward, leading to the generalized hyper- 
geometric function m-\Fn-i- Even for three generations, 2 F 2 is a complicated function. 
From Eq. ( |3.2j ), corresponding to the center of the sun, e~ u is about 3000; for arguments of 
this order of magnitude, it is not practical to calculate the various 2 F 2 of Eqs. (|3.28|) - (|3.30|) 
using the Taylor series expansion. 

Sections IV and V are devoted to the issue of how these 2 F 2 can be calculated. Of the nine 
2 F 2 that appear in Eqs. ( |3.28| )-( p.30| ), three of them are evaluated approximately in Sec. IV. 
In Sec. V, the remaining six 2 F 2 are then expressed in terms of another set of generalized 
hypergeometric functions 3 iq. These 3 iq have convenient integral representations, useful for 
numerical evaluation. The approach described in the rest of this paper is the best we have 
found, but there may well be other methods that are superior. 
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IV. 2 F 2 THAT APPEAR IN V 
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as given by Eqs. (gg§)-(ggg). Define 



In this section, we study the three 2 F 2 that appear in the expressions of \ ipf* 1 and 
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Let 2 -F 2 (l) be considered in some detail; the treatments of 2 -F 2 (2) and 2 -F 2 (3) are similar. 

From the definition of 2 F 2 by its Taylor series, which is always convergent, it is straight- 
forward to verify, by closing the contour of integration in the right half-plane, that 
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This integral is to be evaluated approximately for u large and negative. Except for an 
indentation near s — 0, the contour of integration can be taken along the imaginary axis. It 
is therefore convenient to let 
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when z is purely imaginary. Using ( |4.8| ), the integrand in ( f4.4| ) is 
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2 (t) = -r(ln |r| - 1) - (r + Pl )(ln |r + pi| - 1) + (r + p 2 )(ln |r + p 2 | - 1) 
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Eq. ( |4 . 1 3| ) together with the fact that Q\ (r) is a non-decreasing function of r implies that the 
relevant point of stationary phase occurs when r is positive. Let r be the unique positive 
root of the cubic equation 
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then r is the point of stationary phase that gives the leading contribution to 2 F 2 (1). The 
result is, with all gamma functions replaced by the corresponding Stirling formula (fOJ). 
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These are the desired results. 
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V. 2 F 2 THAT APPEAR IN AND ifjf ) 

The next step is to study the 2 F 2 that appear in the expressions (|3.28|) - (|3.30D for ip[ l \ ip 2 ^ 
and ^3 with i = 1,2. It should be emphasized that, while the treatment in the preceding 
Sec. IV involves approximations, what is to be carried out in this section is exact. Since 
exact manipulations are usually more straightforward, the description will be relatively more 
brief in this section. 

Consider the third-order ordinary differential equation 
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A comparison with Eq. (|3.22|) shows that 
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The three linearly independent solutions of Eq. ( |5.1| ) are 
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where f^ 3 \z) is the function treated in the preceding Sec. IV. 

On the other hand, if in Eq. fl5.1|) the independent variable is changed to 



then the differential equation is 
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In this form, two of the three linearly independent solutions of Eq. ( |5.1|) are (in terms of z) : 



fW(z) = z-°* 3 F 1 



P\z) = z- 



«2 



3^1 



-(3\ + at!, 
1 — a 2 + cei 

-(3! + a 2 , 
1 — «i + a 2 



-(3 2 + «i, -f3 3 + a x 
-(3 2 + a 2 , -f3 3 + a 2 



-i 



(5.10) 
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For a discussion of the generalized hypergeometric functions p F q with p > q + 1, see for 
example Chapter V of 

It should be added parenthetically that, for the present purposes, z _/31 , as an example, 
can sometimes be very small because by Eqs. (|3.21| ) and (|5.2|) 



-Pi 



(5.11) 



Thus, when /ii is large, which is often the case, for numerical calculations it is convenient to 
replace z~^ by e~ lfllU , and similarly for the other powers of z multiplying the 2 F 2 and 3F1. 

Since Eq. ( |5.1|) is a third-order linear ordinary differential equation, the two solutions 
as given by Eq. ( p,10| ) must be expressible as linear combinations of the three solutions of 



Eq. ( |5.7|) . The result is 

f {1 \z) = 



7iT(l - a 2 + ax) 



r(-/3x + ax)Y{-(5 2 + ai)r(-/3 3 + a x ) 
r(-/3 2 + A)r(-/5 3 + A) 



1 



+ 



+ 



T(l — q;i + /9i)r(l — a 2 + (3\) sin 7r(ai — (3\ 

r(-/3 1 + /3 2 )r(-/3 3 + /3 2 ) 1 



r(l - a x + /5 2 )r(l - a 2 + /3 2 ) sin 7r(«i - /3 2 ) 

r(-A + /? 3 )r(-/? 2 + /? 3 ) 1 

T(l - ai + /?3)r(l - a 2 + sin 7r(ai - /3 3 ) 



/ (2) W 
/ (3) (^) 



(5.12) 



ai»->«2 



(5.13) 



Eqs. (|5.12 ) and Q5.13]) can be considered to be the formulas that express f^'iz) and 
f®(z) in terms of f {1) (z),J i - 2 \z) andj^(z). Since f^(z) is given by Eqs. (jOp-([|T7[), 
it only remains to obtain p l \z) and f^(z), which are given by 3-F1 rather than 2 F 2 . 

The advantage of 3-F1 over 2 F 2 is due to the integral representation 



b 



a 2 , a 3 



x 



r (ai 



.7; 



«i 



/•oo 

/ dte- xt t a ^ 1 2 F 1 {a 2: a 3 ;b;-t) : 
Jo 



(5.14) 



where 2 -F\ is the usual hypergeometric function. This Eq. ( |5.14j ) follows from Eq. (8) on 



p. 214 and Eq. (1) on p. 215 of [11|, and has the great advantage of not containing in 
the integrand any gamma function that depends on the variable t, because such gamma 
functions require excessive computer time. 

For the 2 F\ in the integrand of Q5.14 ), the representation 



2 Fi(a,6; c; -t) 



-T(c)e 



4r(6)r(c — b) sin irb sin 7r(c — b) 

) c - h -\l + ts)- a ds 



x 



s b -Hl 



V 



(5.15) 



is convenient, where V is the Pochhammer contour of Fig. 2. In the numerical evaluation 
through Eqs. ( |5.14 ) and ( 5.15|) , it is convenient to deform the contours of integration through 
the point of stationary phase. This point is discussed in more detail in the Appendix. 
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FIG. 2. Pochhammer contour used for the evaluation of the hypergeometric function 2-^1- The 
point of stationary phase, denoted sq, occurs along the path from NE to SW. 



VI. SUMMARY AND DISCUSSION 

The Mikheyev-Smirnov-Wolfenstein effect in the sun has been studied and the differential 
equations solved for three types of neutrinos coupled through their mass matrix. Under the 
assumption, as expressed by Eq. (|2.6|), that the electron density can be approximated by an 
exponential function, the MSW differential equations are solved exactly in Sec. Ill, especially 
Eqs. ( |3.27| )-( |3.3U| ), in terms of the generalized hypergeometric function 2-^2 • The method 
used there can be immediately generalized to N types of neutrinos coupled in the same way, 
leading to at-i-Fat-i- 

Since 2F2 cannot be considered to be a familiar function, Sees. IV and V are devoted to 
the issue of how they can be calculated. Of the nine 2-P2 that appear in Eqs. fl3.28|) -( [3.30D , 



three are evaluated approximately in Sec. IV. In Sec. V, the remaining six 2F2 are then 
expressed in terms of another set of generalized hypergeometric functions, 3F1. These 3F1 
have convenient integral representations given by Eqs. ( |5.14| ) and ( |5.15|) , useful for numerical 
evaluation. Eqs. (|5.14j) and ( |5.15|) can also be used to get asymptotic expressions for 3 F 1 
and hence 2-^2, but this development is not discussed here because it is not needed for the 
study of neutrino oscillations [[[]. 

Finally, it should be mentioned that, while the calculation of Sec. IV can be generalized 
immediately to N types of neutrinos, there does not seem to be a similar straightforward 
generalization for Sec. V. More precisely, while n-iFn-i can be related to nF^-2, the 
generalization of Eq. (|5.14 ) involves n-iFn-2, which does not have an integral representation 



similar to Eq. ( |5.15|) involving a single integral. We believe that, for ^Fn-2, it is necessary 
to use an (N — l)-fold integral to avoid having gamma functions that depend on the variables 
of integration. 
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APPENDIX: 



In this Appendix, Eqs. ( |5.14p and ( p . 1 5[ ) are studied further. For definiteness, let a\ and 
a 2 be defined by (|5.3|), and consider the 3 Fi that appears in the first equation in (|5.10| ). Use 
the inequality (|3.12|) to define the four positive quantities 





= M2 


- ^2, 


V 


= C^2 


- Hi, 


C 


= ^3 


- u 2 , 


c 


= ^3 


-uj 2 <( 



(Al) 



then Eq. ( f5.14j ) gives, for this 3F1 under consideration, 

3-^1 



1 + z£, 1-irj, 1 + z( 
1 + K', 



1 / \ l-fif- Z* 00 



(A2) 



r(i + iO 

where, by Eq. (|3.21|) , 

2/ = e~" (A3) 
is large for w negative and large. The substitution of Eq. (|5.15| ) into Eq. ( |A2| ) then gives 



1 -irj, 1+iC 

1 + iC, 



1 



r(i + iO 



/ i7r/2 \ 



r(i + ic')e 



f\ p -m(l+iC') 



T(l — ir])T(i(r] + £'))4sin7r(l — irj) smiri(ri + (') 



(A4) 



where 



/= [°° dt [ dse-^t^s-^n + sY^-^n + s + st)- 1 -^ 

Jo Jv 

roc r 

= dt ds (1 + s + st)' 1 e 
Jo Jv 



-1 -i4>(t,s) 



with 



(f)(t,s) = ty - £\nt + r}\ns - (( - C) ln (! + s) + CM 1 + s + st). 



(A5) 



(A6) 



This <p is of course not related to the <pj of Eq. ( |2.4| ) . 

For the double integral in Eq. (|A5|) , the point or points of stationary phase are given by 



d<p(t, s) d<p(t, s) 



dt 



ds 



0. 



(A7) 



or 
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V-- + -r-^ = 0, (A8) 

t 1 + s + st ( J 



and 

v C-C , (i + OC 



0. (A9) 
s 1 + s 1 + s + st 

Elimination of s between Eqs. (|A8|) and ( |A9| ) gives the cubic equation for t, 

£(t) = 0, (A10) 

with 

£(t) = y 2 t 3 + y(y + C - 2£ - V )t 2 + [((' - 2£)y - (C - £)(£ + ^ - £(C - 0- (AH) 
When C' — ► C> the three roots of Eq. ( |A10|) are explicitly 

t = t + = (2y)- l {-(y - f - 77) + [(y - £ - r/) 2 + 4^] 1/2 } > 0, (A12) 

t = t. = (2y)- 1 {-(y -Z-v)-i(y-Z-v) 2 + ^y] 1/2 } < 0, (A13) 
* = -(C-0/y<o. (Ai4) 

Since by (jAl|) 

^-) = - T ^ r (C-O>0, (A15) 

we have 

{< 0, for t sufficiently large and negative, 
£;=o; <0 ' < Ai «) 
> 0, for if: sufficiently large and positive. 

It therefore follows that the cubic equation ( |A10| ) has one and only one positive root. Since 
the t integral in Eq. ( |A5|) is along the positive real axis, this implies that there is exactly 
one relevant point of stationary phase, say (t , s ), with t given by this unique positive root 
of Eq. (|A10 ). Furthermore, since 

t(t+) < 0, (A17) 

it follows from ( A16 ) that 

to > t+. (A18) 
From Eq. ( |A8| ) , the corresponding s is given by 

a 



so = - 
and hence 



l + t + 



(A19) 



l + tor 1 < s < 0. (A20) 



For the numerical evaluation of this 3-F1 of Eq. (|A4]), using Eq. ( |A5| ) , it is convenient and 
efficient to choose the Pochhammer contour V to go through this value of Sq, and also to 
deform the contour for the t integration so that at t it is locally the path of steepest descent. 
The other 3-Fi's can be treated in similar ways. 
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